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Abstract 

Several physical systems in condensed matter have been modeled approximating their 
constituent particles as hard objects. The hard spheres model has been indeed one of 
the cornerstones of the computational and theoretical description in condensed matter. 
The next level of description is to consider particles as rigid objects of generic shape, 
which would enrich the possible phenomenology enormously. This kind of modeling will 
prove to be interesting in all those situations in which stcric effects play a relevant role. 
These include biology, soft matter, granular materials and molecular systems. With a 
view to developing a general recipe for event-driven Molecular Dynamics simulations 
of hard rigid bodies, two algorithms for calculating the distance between two convex 
hard rigid bodies and the contact time of two colliding hard rigid bodies solving a non- 
linear set of equations will be described. Building on these two methods, an event-driven 
molecular dynamics algorithm for simulating systems of convex hard rigid bodies will be 
developed and illustrated in details. In order to optimize the collision detection between 
very elongated hard rigid bodies, a novel nearest-neighbor list method based on an ori- 
ented bounding box will be introduced and fully explained. Efficiency and performance 
of the new algorithm proposed will be extensively tested for uniaxial hard ellipsoids 
and superquadrics. Finally applications in various scientific fields will be reported and 
discussed. 
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1. Introduction 

Systems which are composed of many particles can often be modeled as an ensemble 
of hard rigid bodies. Such description is particularly successful when excluded volume 

interactions are dominant and internal vibrational degrees of freedom are negligible. 
Despite the absence of any attraction, particles interacting with only excluded volume 
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interactions exhibit a rich phase diagram with a multiplicity of phases, especially when 
the shape is a non-spherical one[Tl |51 |3J HJ |S] . 

Hard Spheres (HS) are a classical example of hard-body model, which has been 
particularly useful to understand the basic correlation which develops in simple flu- 
ids [H [ll m |9j [To] and provides hints on the slow dynamics which characterize liquids 
approaching the glass transition, where packing effects become even more significant 
[TTJ [T^l US] • Even in the case of molecular fluids, HS models are a good starting point 
for sophisticated liquid matter theories [6J. 

Non-spherical models of rigid bodies are crucial to understand the role of the rota- 
tional degrees of freedom, as well as the role played by the shape in determining the 
system's physical properties. Onsager, introducing a hard sphero-cylinder model (HSC) 
for liquid crystals, showed that, by only changing the aspect ratio of the particles, a ne- 
matic phase can become thermodynamically stable |14j. In Onsager's theory, the internal 
energy of the system is zero and only the entropy, coming from translational and rota- 
tional degrees of freedom of the particles, contributes to the free energy of the system. 
Systems showing an "entropically driven" phase transition have been extensively studied 
over the last 60 years [HI [16]. The study of such transitions has been based on extensions 
of the original Onsager's theory [TTJ [THl [111 [20] > and complemented by experiments (for 
a recent review see \W ) . 

Most of the information for hard-body systems has been calculated using numerical 
simulations [T9j |2l]. Several numerical techniques have been developed in the past to 
simulate particles interacting with only excluded volume interactions. The essence of 
these numerical algorithms involves the evaluation of the overlap between different objects 
or, equivalently, their geometrical distance. The flrst simulations of HS [55] were carried 
out by Alder and Wainwright in 1957, and they provided the first evidence of a crystal 
phase in the case of spherical hard particles (disks and spheres). In 1972 Vieillard-Baron 
|23j published a numerical investigation of the phase diagram of a two-dimensional hard- 
ellipsoids (HEs) fluid, introducing an overlap criterion for HEs suitable for a Monte-Carlo 
(MC) simulation. Building on the work of Viellard-Baron, Frenkel et al. [21] investigated 
in 1984 the phase diagram of a tridimensional system of HEs, through MC simulations. 
Perram and Wertheim [2T introduced a simpler overlap criterion, which has been recently 
used by Donev et al. to perform molecular dynamics simulations of HE [25j [261 [27] and 
which has been also generalized to any couple of smooth convex shapes [28 ] 129 ] [30 ] I31j. 
Simulations of particles with more complex shapes have also been reported. For example 
in 1986 Stroobants et al. carried out MC simulations of a system of hard sphero-cylinders 
(HSC), i.e. molecules consisting of a hard cylindrical rod of length L and diameter D, 
capped at each end by hard hemispheres also of diameter D. These simulations are 
computationally less demanding |32j than the HE ones. More recently, MC simulations 
of hard-cylinders (HCY) have been performed [33] in order to look for a cubatic phase, 
which had been reported for cut hard spheres [33]. The oblate version of HSCs (i.e. 
the solid resulting from the intersection of two spheres) has been investigated by MC 
simulations under the name of UFO (the name comes from the shape of the particles) 

m- 

Molecular dynamics simulations of hard bodies are less common than the Monte 
Carlo ones, since the implementation of the overlap criterion between hard bodies must 
be complemented with an algorithm estimating the collision time between them. The 
evolution of the system in conflguration space is propagated from one collision to the 
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next, giving rise to the so-called event-driven molecular dynamics (EDMD). In EDMD 
the predicted collision times between hard bodies are computed and stored into a time- 
ordered event calendar (as it was first done for HSs by Rapaport although different 
techniques also exist in literature |37j). Such technique has been recently extended to 
Brownian Dynamics of HSs [38l 139] . 

A basic scheme for simulating hard non-spherical bodies is based on the standard 
MD algorithms, where at the end of each time step a check for possible overlaps is 
performed and the simulation is "rewinded" when an overlap occurs |401 141j . This scheme 
is obviously inefficient. An overlap potential of two hard bodies A and B is a function 
F{A, B) , such that F <Q\i A and B overlap. In general it is not easy (and not necessarily 
possible) to find an analytic form for the overlap potential of two rigid bodies of arbitrary 
convex shape, except for the aforementioned cases of UFO, HSC, HCY, HCS and HE. 
Overlap potentials for hard rigid convex bodies can be found in |3]. 

An EDMD algorithm for non-spherical objects employing an event-calendar has been 
recently proposed by Donev et al. and tested on several systems of hard particles 
[SSllMlEniE] ■ Such algorithm [25l|26] requires the use of overlap potentials [24','2y , like 
in MC simulations. In the present paper a different route to provide a general algorithm 
to simulate hard particles will be presented. 

If the surface of two hard rigid bodies (HRB) is smooth enough (i.e. first and sec- 
ond order derivatives are defined over the whole surface) a possible overlap potential 
is provided by the minimum distance between the two surfaces (defining the distance 
as negative when two rigid bodies overlap). In this paper we illustrate a method to 
calculate the distance between two generic convex HRBs based on a Newton-Raphson 
(NR) method. We will also illustrate a simple algorithm to efficiently evaluate a guess 
of the collision contact point and time. Starting from this guess true collision point and 
timecan be calculated by solving a reduced system of equations, again through a Newton- 
Raphson method. It is worth noting that this algorithm can handle grazing collisions, 
i.e. collisions in which two rigid bodies touch slightly [26' , without tuning the algorithm 
parameters significantly. Furthermore, in the case of HEs and superquadrics (SQs), we 
illustrate a new kind of neighbour list based on oriented bounding boxes that can also 
be easily generalized to more complex shapes. Like the algorithm proposed by Donev et 
al. |30} I31j , our method can be easily generalized to arbitrary convex shapes (but also 
decorated with localized patchy interactions, see below) with similar efficiency. 

The present algorithm has been already applied successfully to the simulation and to 
the study of various systems. It has been implemented in the study of structural and dy- 
namical properties of uniaxial HE [HI With this code, adding the possibility of hav- 
ing localized patchy square- well interactions, it has been possible to study the statics and 
the dynamics of primitive models of Water [44^ and Silica [45, , as well as the irreversible 
gelation of a model inspired by stepwise polymerization of bifunctional diglycidyl-ether 
of bisphenol-A with pentafunctional diethylenetriamine [151 H71 H5] . More recently this 
code has also been generalized in order to simulate a coarse-grained model of biological 
systems and primitive models of proteins. 

In Section |2] we introduce the general algorithm for the EDMD of rigid bodies. This 
requires discussing the HRB motion, the evaluation of the distance and the time of 
collision among HRBs and the optimized linked list method. In Section [3] we specialize 
the results of Section |2] to the specific case of HEs; in particular, we discuss a new 
nearest-neighbours list method, that over-performs the simple linked lists method usually 
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employed in EDMD in the case of very elongated HRBs. In Section |4] the performance in 
the case of HEs at various densities, elongations and with respect to the main parameters 
of EDMD is analyzed, while in Section |5] the performance of the algorithm in the case of 
SQs is investigated. In SectionlGlsome perspectives and applications of this new algorithm 
are discussed, and in Section [7|conclusions are drawn. 

2. An event-driven algorithm for rigid bodies. 

2.1. Geometry of rigid bodies. 

The orientation of a HRB can be represented by the 3 column eigenvectors (with 
i = 1,2,3) of the inertia tensor expressed in the laboratory reference system. These 
vectors form an orthogonal set and can be arranged in a matrix R, i.e. 

R=(uoUiU2)^ (1) 

where indicates the transpose of A. 

This matrix will be referred to as the "orientational matrix" in the following. The 
orientational matrix relates the coordinates x in the laboratory reference system to the 
coordinates x' in the HRB reference system via: 

x' = Rx (2) 

The following discussion focuses on HRBs with finite volume and bounded surface. 
We assume that the equation of the surface of the HRB, in the "HRB reference system" 
with origin in the center of mass and axes parallel to the vectors {u^ji, is of the form 
/(r) — 0, where / changes sign passing from the interior to the exterior of the HRB. 
Moreover, the normal ^ and its first order derivatives are assumed to be properly 
defined. 

2. 2. Motion of rigid bodies 

In the following equations we assume that the three eigenvalues of the inertia tensor 
are all equal to /. The formulas for the free rotation of a symmetric-top case are just 
slightly more elaborated [39] , although the free rotation for a general rigid body involves 
the calculation of special functions @9l [50] and requires a more sophisticated algorithm 
which has been implemented only recently [HI [52]. Nevertheless, this paper being fo- 
cused on an algorithm which predicts collision between HRBs, for the sake of simplicity 
the discussion will be restricted to the case of a fully symmetric inertia tensor. It is 
straightforward to generalize the approach to a completely asymmetric inertia tensor. 

From the angular velocity w = (wxTWy^Wz) of a free rigid body the antisymmetric 
matrix can be built, i.e.: 

(0 -Wz Wy \ 

Wz -w., (3) 
-Wy J 

It is possible to relate the orientation R(<) to the orientation at time t = via [^I53| : 



R(i) = R(0)(I + M) 
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(4) 





(b) 



Figure 1: Examples of solution of Eq. [9] (a) Distance between tow HRBs with a possible spurious 
solution, (b) For non-convex objects, there can be multiple solutions. 



where M is the following matrix: 



1 -C0SM) ^2 



(5) 



with u) = ||w||. Note that if w = then R{t) = R(0). 

The update of position and orientation of a free rigid body is therefore accomplished 

by: 

x(i) = x(0) + vt (6a) 



R(t) = R(0)(I + M) 



(6b) 



where x(i) is the position of the center of mass of the rigid body at time t and v is its 
velocity. 



2.3. Distance between two rigid bodies 

The present algorithm is based on a calculation of the distance between HRBs. In 
the following we will show how such a distance is calculated. Consider two rigid bodies, 
A and B, whose surfaces are implicitly defined by the following equations: 



/(x) - 
5(x) = 



(7a) 
(7b) 



It is also assumed that if a point x is inside the rigid body A then /(x) < ((/(x) < 0), 
while if it is outside /(x) > (.g(x) > 0). The distance d between these two objects can 
be defined as follows: 

c?= min ||xa-xb|| (8) 

The latter equation means that the quantity \\xa — xb|| has to be minimized with the 
constraints that points x^ and x^ belong respectively to the surfaces of A and B. Hence, 
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introducing the Lagrange multipliers a and /3, the distance between the two HRBs A 
and B can be defined as the solution of the following set of 8 equations: 



-a^3xB 



(9a) 



/(xa) = 
3(xb) - 



xb 



where = {x a-,V A-, z a) , xb = {xB,yB,ZB), /xa = ( 



(9b) 
(9c) 
(9d) 



Eqs. (9b I and (9c I guarantee that and x^ are points on A and i?, Eq.(9a 



guarantees that the normals to the surfaces are anti-parallel, and Eq.(9dl guarantees 



that the displacement of ha from x^ is coUinear to the normals to the surfaces. 

Equations ^ define extremal points of d [51] ; for example for two convex HRBs local 
maxima can be found (see Fig. [ija)) and for two general non-overlapping non-convex 
HRBs these equations can have multiple solutions (see Fig. [ijb)), although only the 
smallest one is the actual distance. Therefore, to solve these equations iteratively it is 
necessary to start from a good initial guess of (x^i, xs, a, /3) to avoid finding spurious 
solutions. 

In addition we note that if two rigid bodies overlap slightly (i.e. the overlapped 
volume is small) there is a solution with /3 < 0, that is a measure of the inter-penetration 
of the two rigid bodies; such a solution will be referred to as the "negative distance" 
solution (Fig(2]). 

Finally, we define the quantities di as follows: 

d,^||x^-x^|l (10) 

where (x^, x^, a^, (3i) is a solution of Eqs. the distance d between two HRBs can be 
formally written as: 

d = sgn(/3„i„) min{di} (11) 

where sgn(a;) is the sign function and (3rnin is the (3i corresponding to the solution with 
the smallest d, . 



2.3.1. The Newton-Raphson method for calculating the distance 

The set of equations (|9| can be conveniently solved by a Newton-Raphson (NR) 
method This method, as long as a good initial guess is provided, reaches the 

solution very quickly thanks to its quadratic convergence |55j . If one defines: 



/ /XA + 



F(y) = 



/(xa) 
3(xb) 
\xA + /3/xA -xs/ 



(12) 



^Note that the gradient is intended to be a row vector. 
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d >0 



d<0 




(a) 




Figure 2: When convex objects slightly overlap, the solution of Eqs.(j9| changes sign. 



The Eqs. (|9| become: 

F(y) = (13) 

where y = (x^, xb, a, /3). 

Given an initial point yo, a sequence of points converging to the solution can be built 
as follows: 

ym=y,:+J"'F(y,) (14) 
where J is the Jacobian of F, i.e.: 



/ 




0^ 


2a5xB 



o\ 








T 
.Kb 








V 




I 








(15) 



with I being the identity 3x3 matrix and the null column vector. 

The NR method may not converge if the initial guess is not close enough to the 
root, hence in general it may be convenient for the sake of numerical robustness to use a 
globally convergent NR method (see [SS]) that converges to the solution from any starting 
point. An alternative route to ensure the appropriate robustness is to provide an accurate 
initial guess, e.g. making use of a steepest descent method, as it will be shown in the next 
section. The matrix inversion, required to evaluate J~^, can be performed by standard 
LU decomposition j55j. LU decomposition is of order iV^/3, where TV is the number of 
equations (8 in the present case). In Sec. 2.3.3 it will be shown that the above set of 



equations can be reduced to 5, thus reducing the time to invert the matrix by a factor 
-4. 



2.3.2. Initial Guess for the distance: the steepest- descent method. 

As initial guess of the NR it may be convenient to choose the closest pair of points 
on the intersection of the surfaces with the line joining the center of mass of the two 
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(a) (b) 

Figure 3: Initial guess for the distance to be used as a starting point for NR method, (a) Simple initial 
guess for the distance: initial points and xb for the NR are obtained considering the interceptions 
of the line joining the two HEs centers with their surfaces, (b) A better guess for the distance in the 
case of moderately elongated HRBs (see discussion in Sec. |3.2|for the HEs case). 



HRBs (Fig.|3^). Such guess is reasonable if the objects are almost spherical; otherwhise 
a refinement is needed. 

A general method to refine such initial guess is to minimize the function 

^(y) =«sd||xa-xb||2 (16) 

where y = (x^, x^) with the constraints /(x^) = and (/(x^) = 0, and asD is parameter 
that can be tuned to optimize the steepest-descent (SD). For solving such a problem, 
steepest descent steps are used, followed by corrections that hold the points and xs 
on the surface of the two bodies. The algorithm is the following one: 

1. Choose an initial guess for x^^ and x^ 

2. Evaluate the gradient Dy of D{y): 

Dy = (h^, hs) = 2asD{xA - x^, -(x^ - x^)) (17) 

3. Calculate the components of and hs parallel to the surface, i.e.: 

4 = (hi, hi) (18) 

where 

h|i = V - ■ n^)"/^ (19) 
with ^ e {A,B} and = /x^/||/x^||, hs = gxa/Ugx^ll 

4. Move the two points in the direction of Dy-. 

y' = (x:4,xij)=y-A5i54 (20) 

whith < XsD < 1- 

8 



5. Add a small displacement dy = (Ca/xaj^-B5xb) to the vector y', i.e. 



y" = (x:;,x^)^y' + dy (21) 

where and S^b are such that the constraints are satisfied and the two points 
and yJg belong to the surfaces of the two rigid bodies, i.e. 

/(x:4 + e^/x J - (22a) 

g(x'B + es5x J = (22b) 

6. Terminate if the angle between n and Dy is small enough, otherwise go back to step 
1. 

This procedure provides a guess for -ka and x^. Note that the adjustment of the position 
of the points A and B to hold them respectively onto the surfaces / and g can be 
implemented again with a one-dimensional Newton-Raphson method, that will generally 
converge in few steps if the points are not too far from the surfaces. 

The NR method for the distance requires also a guess for a and (3 (introduced in 
Eq. Q); a = H/xAll/llffxisll, /? = proved to be a good guess. If the accuracy required 
for the convergence of the SD is high enough, the NR method will always converge to 
the correct solution. However, being the SD method much slower than NR, a trade-off 
between accuracy and speed is needed. 

2.3.3. Reduced system of equations 

The system in Eq. ^ can be reduced to 5 equations eliminating Eq. ( 9d ) : 



/x^ + a'5x« (XA + /3/x^ ) = (23a) 

/(xa) = (23b) 
<7(x^ + /3/x^) = (23c) 

In this case the Jacobian is: 

J - fL (24) 

where A is the following matrix: 

A^|^ + /3?^|^ (25) 
oxb oxb ox a 

Evaluation of J^^ can be done again by the LU decomposition. 
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2.3.4- Relaxed Newton-Raphson method 



Eq. (14 1 has been derived using a first order expansion of F{y), but if the matrix J 



is singular or nearly-singular in the proximity of the solution, even with a good initial 
guess the NR step can become too large and the NR method unstable. To solve this 



problem, a damping factor into Eq. ( 14 1 can be conveniently introduced: 



(26) 



where < < 1 and a suitable has to be found in order to stabilize the NR method. 
First, consider either the system of equations ^ or (23 1; let dyi — (dx^, dxB,da'^,dP) 
be the step predicted by NR. The relative variations of a^, f3 can be bounded, imposing 
the conditions ^i|c?a^| = ^jjadaj = e^a^ and ^i|c?/3| = er\P\ with < e,. < 1. 

Taking the modulus of both sides of Eq.(9a) one gets = || /x 4 11/11.9x^5 II > obtaining 
the condition: 

er-ll/x^ll 



6 = 



2||5xJ||a| 



(27) 



Similarly, from Eq. (9d), one obtains the equation ||x^ ^x^H = H/xaIII/^I- Assuming 
that \\xA — xsll ^ I, where I is a typical distance of the system, provides the condition: 



Therefore, the damping factor can be chosen such that: 

^r^ \ \f:x.A I 



^^"""M/3|ll/xJI'2||gx«, 
to ensure that the variations of a and f3 remain sufficiently small. 



(28) 



(29) 



2.4. Prediction of the collision time 

In an event-driven (ED) algorithm one needs to predict the time of collision of pairs 
of HRBs. This means to evaluate — given two objects at time t = and the distance as 
a function of time d{t) — the smallest time tc > such that d{tc) — 0. If a good guess 
of the contact point and time is provided, solving a proper set of equations through a 
NR method (as it will be shown shortly) allows to find the exact contact point and time. 
In the following it will also be shown how to calculate in a very efficient and simple way 
such an initial guess, exploiting the evaluation of the distance between two rigid bodies. 

2.4.1. Newton-Raphson method for the contact time 

The point and time of collision satisfy the following equations: 

/x(x,t) + a2^x(x,t) = (30a) 

/(x,t) = (30b) 

g(x,t) = (30c) 

where / and g now depend also on time because the objects move. Again, it is appropriate 
to employ the NR method to solve such a system; the Jacobian for the NR in this case 
turns out to be: 
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fl /, (31) 

\ gl fft / 

where ft = ^ and <7t = §f • As discussed in [56], such method becomes very unstable 
even for simple convex objects unless one starts from a very good initial guess. To 
construct such a guess, one can find a good bracketing of the collision time, i.e. two 
times ti, t2 such that: 

1. d{ti) > and d{t2) < 

2. <2 ~ ti is small enough to have at the most one collision within {ti,t2) 

3. d{ti), d{t2) are small enough in order to avoid any instability of the NR method (see 
above) . 

Once the latter bracketing has been found, the initial time for the NR is set to ti, while 
a good initial guess for the contact point will be halfway on the distance between the 
two bodies at ii. 

2.4-2. Bracketing the contact time 

To first bracket the contact time a "bounding sphere" (BS) may be used (FigQ. The 
BS of a given HRB A with center is the smallest sphere centered at that encloses 
A. Given two rigid bodies A and B at time t — and their BSs Ca and Cb, one must 
indeed consider three possible cases: 

1. Ca and Cb do not overlap and they will not collide: in this case A and B won't 
collide either. 

2. Ca and Cs do not overlap but they will collide: in this case one has to search for 
a collision within the time interval [^1,^2] where ti is the time when the two BSs 
collide and start overlapping and t2, which is greater than ti, is the time when they 
just cease overlapping. 

3. Ca and Cb overlap: in this case one has to search for a collision within the time 
interval [ti,i2] where ti = and t2 > ti is the time at which the two BSs stop 
overlapping. 

Since the collision prediction between hard spheres (i.e. BSs) is extremely fast this 
technique improves the performance by reducing the number of collisions to check and 
provides a first bracketing (ti, t2) of the collision time. 

To further improve this bracketing, the following overestimate of the rate of variation 
of the distance can be established: 

d{t) < ||vA - vbII + \\^a\\La + ||wb||Lb (32) 

where the dot indicates the derivation with respect to time; wa and are the velocities 
of the centers of mass of A and B\ and are the angular velocities of A and B, 
and the lenghts La, Lb are such that 

La> max |||r' — r 4 III (33a) 

/(r')<0 
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Figure 4: Minimal bounding sphere for a rigid body. 



Lb > max {||r' - r_B||} 

9(r')<0 



(33b) 



the case d(t) > 0. In this case, for a positive time increment dt > one has: 

\liAit + dt) -XB{t + dt)\\ - ||xA(i) -XB(t)|| 



where r^, are the centers of mass of the two rigid bodies. To derive Eq.(|32|, consider 

(34) 



m\ 



dt 



where x^ and x^ are two points belonging to the two rigid bodies such that: 

d{t) = \\^A{t)-^Bm (35) 

If Xyi(i) and x^(t) are the velocities at time t of the points x^(t) and XB(i) respectively, 
it results that: 



||x^(i + dt) - ^B{t + dt)\\ < ||x^(t) - XB(t)|| + \\i±Ait) - ±Bit))dt\\ 



(36) 

because of Eq. ([8|, that defines the distance. From Eq. ( [36| and Eq.([34| the following 
inequality is obtained: 

\d{t)\ <\Wa- VijII + ||wA X (x^(t) - rA{t))\\ + ||wB X (xB(t) - rB(t))|| (37) 

and Eq. p2| results from the fact that: 

\\^A{t)-rAit)\\<LA (38a) 

\\^B{t)-rB{t)\\<LB (38b) 

The case d{t) < is similar. Indeed, performing the substitutions t + dt ^ t' and 
dt —dt' into Eq. (34 1, Eq. (32 1 is obtained again. 

With such an overestimate of d{t) , that from now on will be called d^ax i an efficient 
strategy for the refinement of the bracketing (ii, ^2) of the collision time and point is the 
following one: 
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3 



Set t = ti. 

Evaluate the distance d{t) at time t. 
Choose a time increment At as follows: 




if d{t) > €d 



otherwise. 



(39) 



where ^ inm{LA, Lb}- 

4. Evaluate the distance at time t + At. 

5. If d{t + At) < and d{t) > 0, then ti = t and t2 = t + At, find the collision time 
and point via NR and terminate. An initial guess for this NR may be obtained 
through a quadratic interpolation of the points {t, d{t)), {t + At/2, d{t + At/2) and 
{t + At,d{t + At). 

6. If both < d{t+At) < and < d{t) < ej,, a, "grazing" collision may occur between 
t and t + At (because distance may be first diminishing, then growing). To look 
for a possible collision, evaluate the distance d{t + At/2) and perform a quadratic 
interpolation of the 3 points {{t, d{t)), {t + At/2, d{t + At/2), [t + At, d{t + At)). 
There are 3 possible cases: 

(a) The parabola has a minimum and d{t„i) < 0, then set ti — t and ^2 = tm, find 
the first zero tz of this parabola and use tz and d{tz) to start a NR to find the 
collision time and point, and terminate. 

(b) If the resulting parabola has a minimum between t and t+At, but d(tm) > 0, 
then find the minimum tmb of d(t) between t and t + At with the maximum 
possible accuracy (using a 1-dimensional Brent's method for finding minima). If 
d{tmb) < then set ti = t and t2 = tmb and find the first zero tz of d{t) between 
ti and t2 (with a moderate accuracy, because tz will serve only as initial guess 
for the NR). Use tz and d{tz) to start a NR to find the collision time and point 
and terminate. 

(c) If the parabola does not have a minimum then keep searching (i.e. go to step]?]). 

7. Increment time by At, i.e. t t + At. 

8. if t > t2 terminate (no collision has been found). 

9. Go to step 2. 

Note that if starting at time t the two rigid bodies collide at time tc > t and d{t) > e^, 
then our over-estimating procedure enforces At < tc — t. If the chosen is small enough 
the distance d{t) has only one minimum within the time interval [t, t + At] and the 
above scheme for predicting HRBs collisions ensures that "grazin" collisions are properly 
handled within machine accuracy, i.e. all "grazing" collisions are correctly predicted. 
According to step 6) in the above scheme, indeed, if d{t) > and d{t + At) > 0, a 
quadratic interpolation is used to search for a negative minimum (i.e. a minimum such 
that d{t„i) < 0). If such a negative minimum can not be found, an attempt to find it 
is made with the maximum possible accuracy using a suitable numerical method (e.g. 
Brent 's method) . We stress that must not be chosen with unreasonably tight tolerance 
not to miss grazing collision within machine accuracy and, although finding the minimum 
with Brent's method can be time-consuming, this method is only a second option after 
the quadratic interpolation, which is faster and finds the minimum in the majority of 
cases. 
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2. 5. Elastic collision of two hard rigid bodies 

As two rigid bodies collide, one has to evaluate the new velocities of centers of mass 
and the new angular velocities after the collision. Imposing the conservation of impulse, 
angular momentum and energy, the velocities after the collision can be evaluated as 
follows: 

-> + TO;4^ApAsn (40a) 

vs ^ vs - rrig^ApABn (40b) 

WA ^ V/A + ApAslA^irA- xc) ^ n (40c) 

ws ^ w_B - ApABlB^i^B - xc) X n (40d) 

where is the contact point, n is a unit vector perpendicular to both surfaces at Xj^, 
I A, Ib are the moments of inertia of the two colliding rigid bodies, m^, are their 
masses and 

2(vA + WA X (xc - r^) - vs - ws X (xc - rs)) 

lApAB — 31 zr\ 31 (41) 

+mg +/^ ||(r^ -xc) X n|| +/g ||(rB - xc) x n|| 

2.6. Linked lists for bounding spheres 

Predicting collisions is the most computationally demanding part of an EDMD. In 
order to speed up a EDMD of hard spheres, one can use linked lists [35] to avoid to check 
all the N"^ possible collisions among N objects. In the linked list method, the simulation 
box is partitioned into cells and only collisions between particles inside the same 
cell or its 26 adjacent cells are accounted for. This means also that, whenever an object 
crosses a cell boundary going from cell a to a new cell 6, it has to be removed from cell 
a and added to cell b. 

As a first step, one can recover the same method using the BSs location as particles. 
In this case, the cubic box of side L containing N identical rigid bodies is divided into 
cells so that each cell has side length of the order of the BS diameter ctc. After that 
linked lists of these BSs are built and updated as in a ordinary EDMD of hard spheres 
[55] . To predict collisions of a rigid body A, one takes into account only the rigid bodies 
that have their BSs inside the 26 adjacent cells or in the same cell as A (see P6 for more 
details) . 

Unfortunately, BSs are not very efficient if the volume of the BS is much bigger than 
the volume of the rigid body, as it happens for example for ellipsoids of high elonga- 
tion [25] [55]. In this case, the number of possible collisions which must be calculated 
makes such procedure computationally inefficient. In Section [3. 3 j an alternative and more 
efficient method for objects with high elongations will be illustrated. 

2. 7. Putting all together: hard rigid bodies event-driven algorithm 

In an ED algorithm many events may occur, such as collisions between particles, 
cell crossing (if one uses linked lists), saving of a system snapshot, output of measured 
quantities, etc. All these events should be ordered in a calendar so that the next event 
to happen can be easily retrieved; at the same time, insertion and deletion of events in 
the calendar should be done as quickly as possible. 

One elegant approach was introduced almost thirty years ago by Rapaport [361 , who 
proposed to arrange all the events into an ordered binary tree (calendar of events), so that 
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insertion, deletion and retrieving of events could be done with an efficiency O^logAf), 
0(1) and 0{logN') respectively, where Af is the number of events in the calendar. This 
solution is adopted to handle events calendar in our simulations; all the details of this 
method can be found in [36 . 

Considering that all the tools to develop a standard ED algorithm for hard rigid 
bodies have been illustrated, the algorithm can be resumed as follows: 

1. Initialize the events calendar (predict collisions, cell-crossings, etc.). 

2. Retrieve next event £ and set the simulation time to the time of this event. 

3. If final time has been reached, terminate. 

4. If f is a collision between particles A and B then: 



(a) change angular and center-of-mass velocities of A and B according to Eqs. (40) 

(b) remove from calendar all events (collisions, cell-crossings) in which A and B are 
involved. 

(c) predict and schedule all possible collisions for A and B. 

(d) predict and schedule the two cell-crossings events for A and B. 

5. If £ is a cell crossing of a certain rigid body A: 

(a) update linked lists accordingly. 

(b) remove from calendar all events (collisions, cell-crossings) in which A is involved. 

(c) predict and schedule all possible collisions for A using the updated linked lists. 

6. go to step 2. 

The novel aspect of the present algorithm is in the time-consuming step (4)c, for 
which we have provided in the previous section the implementation details. 



3. Hard ellipsoids with an axis of symmetry 

In this Section the details about a specific case for which the algorithm has been 
tested will be provided [12] : hard ellipsoids with an axis of symmetry. Such ellipsoids are 
characterized by the elongation Xq, i.e. the ratio of the length of the symmetry axis with 
respect to any of the other axes. A new efficient implementation of nearest neighbor lists 
(NNL) for hard ellipsoids will be also discussed and a careful test of the performance of 
this approach will be shown. 

3.1. Evaluation of the distance between two ellipsoids 

The surface of a hard ellipsoid can be implicitly defined as follows: 

(x - r^)^X^(x - r^) = (42) 

where r^ is the position of the center of mass of the ellipsoid and is a 3 x 3 positive 
definite matrix. In particular if at time t = the rigid body reference system coincides to 
the laboratory reference system it turns out that the free evolving hard ellipsoid surface 
is: 

Xa - Rlit)XA{0)RAit) (43) 

where 




and a, b, c are the three semi-axes of the ellipsoids. For hard ellipsoids with an axis of 
symmetry, the values of two semiaxes are equal. 

To evaluate the distance between two ellipsoids A and _B at a time t one has to 
solve either Eq. ^ or Eq. (23) by the Newton-Raphson method. For Eqs. ^ in this 
particular case the Jacobian becomes: 
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(45) 



where = 2Xa(xa - r^) , ub = 2Xb(xb - r^) and 



X^ = R^(t)X^(0)R^(t) 



(46) 



with ^ e {A,B}. 



For Eqs. (23) the Jacobian turns out to be 

''2X^ 



a^A 2anB 




(47) 



where yia, and X^ are the same as before, xb = x^ + /3n^ and 

c — 2a^ Xb • 

d^2l3nl-XA 
A^2a^XB-{2pXA + I) 



(48a) 

(48b) 
(48c) 



3.2. A better guess of the distance 

As already discussed, the NR method needs a good guess of the starting point and for 
this purpose a steepest descent method has been presented in Sec. |2.3.2[ For ellipsoids 
with moderate elongations (0.2 < Xq < 5.0), the initial guess can be calculated in an 
alternative and very simple way. The simplest possibility is to use as a first guess for x^ 
and xb the intersections of the vector tab — J^a — J^b, joining the centers of mass of the 
two ellipsoids, with their surfaces. This guess is quite rough and a possible improvement 
can be achieved, as explained in the following. 

First of all the components of tab in the reference systems of the two ellipsoids are 
calculated, i.e. 

v'a = Ra ■ vab (49a) 

V^j = Rb • TAB (49b) 

then these two vectors will be scaled by the semi- axes of the ellipsoids and their compo- 
nents will be transformed back to the laboratory reference system, i.e.: 

= RliS ■ ^,'^) (50a) 



VB = Ri;(S-v^) 
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(50b) 



where and are the orientational matrices of A and B respectively and 



S = , I b I (51) 

Finally the intersections with the HE surfaces of the vector centered at the center of 
mass of A and centered at the center of mass of B are computed and these two points 
are the desired guesses for and for the NR. The effects of such trasformations are 
shown in the right panel of Fig. ([3| . 

At this point, distance evaluation can be started with the initial values for a and /3 
set to 0. This method provides a speed-up of around 10%. 

3.3. Nearest Neighbour Lists 

Increasing the elongation of the ellipsoids, the linked list method illustrated before 
becomes progressively less efficient at moderate and high densities [25^ . Indeed given one 
ellipsoid A, when using the linked lists one has to predict the collision times of A with 
all ellipsoids in the cell of A and with all ellipsoid in the 26 adjacent cells. In the case of 
rotationally symmetric ellipsoids, the number of collision time predictions grows as Xq if 



Xq > 1 and as I/Xq if Xq < 1 [26] (see also Sec. 4.3 1. In order to reduce the number of 
predictions at very high or very small elongations, Donev et al. [25J suggest to surround 
each particle A with a bounding neighborhood having the same shape as A (i.e. in the 
case of HEs they use ellipsoids) and to predict collisions only between particles having 
overlapping bounding neighborhoods. Similarly to what is proposed by Donev et al. j25j 
we suggest to build an oriented bounding parallelepiped (OBP), instead of an ellipsoid, 
at a certain time t around each HE and to predict collisions only between HEs having 
overlapping parallelepipeds (Fig.[5|. 

Each parallellepiped encloses the corresponding ellipsoid completely, and it is centered 
at the center of mass of the ellispoid itself. More precisely, given an ellipsoid A with semi- 
axes a, 6, c and center of mass r^, the vertices Vq of the parallelepiped, with a G {1 ... 6} 
are: 

Va = + 0-2(a)(s<xi(Q) + ANNL)'a.a^(a) (52) 

where 

cri(a) = 1 + (a -1)^2 

0-2(0;) = 2 (a mod 2) -1 (53) 

and s ~ (51,^2,53) — (a,5, c). Moreover {ui},jg{i 2,3} are the principal axes of the 
ellipsoid, i.e. they are such that: 

XaUc = SaUa (54) 

and Aatjvl is a positive parameter that can be tuned to optimize the performance of the 



NNL (see Sec. 4.51 



Given an ellipsoid A, the set of ellipsoids having their parallelepipeds overlapping 
with the parallelepiped enclosing A, is the NNL of A. Each parallelepiped i is immobile 
and contains only the i-th ellipsoid; if ti is the time when the ellipsoid will collide with its 
containing parallelepiped, its NNL will have to be rebuilt not after the time tr = miiii{ti}. 
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Figure 5: Nearest neighbour lists are built through OBPs: ellipsoids A and B are nearest neighbors, 
and they may collide among them before colliding with their OBP; neither A or B can collide with C 
during this time. The "thickness" of the OBP is Ajvjvl- 



In addition, if a collision between two ellipsoids i and j occurs, then the new contact 
times ti and tj of these two ellipsoids with their parallelepipeds will have to be evaluated 
and the new time for rebuilding the NNL lists, i.e. t'^. = mm{tr,ti,tj}, will have be set. 



3.3.1. Distance between a rigid body and a plane 

For predicting the time of collision between an ellipsoid and a parallelepiped the 
distance between an ellipsoid and each of the 6 faces of the parallelepiped must be 
calculated. This means that one has to evaluate the distance between the surfaces defined 



through Eqs. (42) and a plane defined through the following equation: 



np • (x - rp) 







(55) 



where both rp and np do not depend on time because the plane is immobile (Fig. |6|. 

Again a NR method can be used to calculate this distance. In Sec. 2 the second 
rigid body B defined by Eq. (7b I may also be a plane (as defined through Eq. (55l) 



or, alternatively, a plane can be thought as a limiting case of a very large ellipsoids, in 
which case the Jacobian needed to solve Eqs. ([9| turns out to be: 
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(56) 



while the Jacobian for solving Eqs. (23 1 is: 





(57) 




Figure 6: The initial NR guess to calculate the distance between a plane P and an ellipsoid E uses the 
points x^; and xp. Here np is the normal to the plane; the center-of-mass of the ellipsoid is on the line 
from Xp with direction np. 



A good guess is required also to solve Eqs. ^ or Eqs. (23 1 through a NR method. Let 
I be the axis passing through the center of the ellipsoid and parallel to the vector np. 
If intersections points and xp of I with the ellipsoid and the plane respectively are 
evaluated, then these two points (if Eqs. ^ are used) or only x^ (if Eqs. (23) are used) 
proved to be a good initial guess for the NR method. 

3.3.2. Prediction of the collision time 

Given a good bracketing of collision time of an ellipsoid A with one of the six planes, 



Eqs. (30), where / defines the surface of an ellipsoid and g the surface of a plane, can 



be solved numerically. The Jacobian to use in the NR is the following: 

^2Xa -2anp 
J = ( k\ (58) 



where h = f2 x — X^iv^i, k = —v^Xyix + x"^ f2 x — x X^v^i and f2 = rj^X^i — Xy^n^i, 
X = X — r^,. Here Q,a is the matrix Q, defined in Eq. ([3]) corresponding to the angular 
velocity of the ellipsoid A and is the center of mass velocity of A. 

As for the collision of two ellipsoids, the collision time can be bracketed using an 
overestimate of d{t) (being d{t) the distance between an ellipsoid and a plane of the 
parallelepiped). If the ellipsoid surface is the locus of points such that /(x) — and the 
plane is defined as in Eq. ( |55[ ), the distance between an ellipsoid and a plane can be 
defined as follows: 

d— min Ixyi • npl (59) 

/(xa)=o' 

where fip = np/||np||. Similarly to what has been done for two ellipsoids, it can be 
proved that 

rf(i) < |vA-np| + LA||wA|| (60) 
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From now on, the latter overestimate will be called dpm- It is now possible to illustrate 
the algorithm to find the collision time of an ellipsoid with its OBP. Consider an ellipsoid 
A and the six planes of the corresponding parallelepiped labelled by {pi}i<i<6- With 
respect to each plane there are six different overestimations of d{t), labeled as <fp^, and 
six different distances di{t) at time t. The algorithm is the following one: 

1. set t = ti. 

2. Evaluate the distances di{t) at time t. 

3. Choose a time increment At as follows: 

(mmi{di{t)/di^}, iidiit)>eT'; 

At= { e!?"' , . (61) 

— 7) otherwise. 

where e^"' < La- 

4. Evaluate the distance at time t + At. 

5. If for at least one i one has di{t+At) < and di{t) > {). then the solution is bracketed. 
Find the collision time and the collision point through the NR. For the initial guess 
of the collision time, use the first zero of the parabola resulting from a quadratic 
interpolation of points {t,d{t)), {t + At/2,d{t + At/2)) and {t + At,d{t + At)). After 

that, terminate. 

6. For all distances such that < di{t + At) < ej, and < di{t) < ed, evaluate the 
distance di{t + At/2), perform a quadratic interpolation of these 3 points ( {t, di{t)), 
{t + At/2, d,{t + At/2), [t + At, di{t + At)) and find if the resulting parabola has 
any zeros. If this is the case, refine the position of the smallest zero using again the 
NR. 

7. Increment time by At, i.e. t — > i + At. 

8. if time is greater than t2, terminate (no collisions have been found). 

9. Go to step 2. 

Grazing collisions in this specific case are dealt with differently from collisions between 

two HEs. The trick is to consider, for predicting the escape time of a HE, a smaller 
bounding box. In particular the chosen distance between each of the 6 planes forming 
the bounding box and the embedded HE is A jv jvl — e""' • With this choice of the bounding 
box dimensions, if a collision is missed, due to a grazing collision, the HE is not outside 
its bounding box anyway and values for e^^^ around 10~^ b can be safely chosen. 

The times ti and t2 are evaluated making use of the BS of A as illustrated in the 
following. Consider the ellipsoid A and its BS Ca at a certain time t when the NNL has 
to be rebuilt. Two possible different cases may occur: 

1. Ca is enclosed into the parallelepiped: in this case ti is the time when the BS will 
collide with one of the six planes of the parallelepiped and t2 > ti is the time when 
the BS stops intersecting the plane. 

2. Ca intersects the parallelepiped: in this case t\ = t and t2 > t\ is the time at which 
the two BSs cease to overlap. 
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3. 3. 3. Overlap of Parallelepipeds 

For building NNL all the parallelepipeds overlapping with a given one have to be 
found. This task can be accomplished quite efficiently, using a technique well known in 
computer graphics |57j . Consider two parallelepipeds corresponding to ellipsoids A and 
B with centers ta and r^, principal axes vl^ and and semi-axes a,b,c. 

Consider the following straight lines originating from the center of A: 

tj ^rA+ (62) 

with 

wj e {{u^}a, {uf }„, {u^ X uf (63) 

where a, P — 1,2,3, ^ £ and j = 1 ... 15 labels all the possible directions. 

The centers of the ellipsoids ta and will be projected onto all these lines, obtaining 
the points f ^ and f ^ and the vertices of the two parallelepipeds will be projected as well, 
obtaining the points pfj and pfj, where i = 1 . . . 8 labels all the possible vertices of a 
parallelepiped. 

Using the projected vertices one can build for each direction j two spheres with centers 
and and radii: 

Pf = niax{||pf^- -f/||} 

pf = max{||p5-ff||} (64) 

The two parallelepipeds do not overlap if and only if all these pairs of spheres do not 
overlap. 



3.4- Prediction of the collision time between two ellipsoids 

Exploiting linked lists and BSs, one obtains a time interval [^1,^2] that brackets 
the possible contact time between two ellipsoids A and B. Making also use of NNL, a 
better bracketing of the collision time can be estimated, with the bracketing time interval 
becoming [^1,^2], with: 

i2 = mm{t2,tr} (65) 

where is the time at which the NNL rebuild is scheduled. Then starting from t = ti 
the collision time can be finely bracketed by applying the algorithm illustrated in Sec. 
for the general case of hard rigid bodies with the substitution i2 ^ ^2 • 
The bracketing of the contact point and time for the collision between two HEs can 
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be used in the NR to solve Eqs. (|30| with the following Jacobian: 

(66) 



^2{Xa + a^Xs) 




where 

h = - x^v^ + a^(ris x_B - xbvb) (67) 

and 



kA = -V^X^XA + XA - XA Xava 
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= ~v|XbX^ + Xb - xb XbVb (68) 

Here, one has: flA = ^a^a — ^a^a, = ^b^b — ^b^b, ^a = x — r^, = x — r^, 
where flA is the matrix ft defined in Eq. ([s]) corresponding to the angular velocity wa 
of the ellipsoid A, and is the same matrix for the ellipsoid B. 



3.5. Event-driven molecular dynamics of Hard Ellipsoids 

The ED molecular dynamics of HEs is similar to what illustrated in Sec. |2.7[ The 
only addition is the use of NNL to speed the simulation up. In the following the scheme 
of an EDMD of HEs is given: 

1. Initialize the events calendar (predict collisions, cells crossings, etc.). 

2. Retrieve next event £ and set the simulation time to the time of this event. 

3. If final time has been reached, terminate. 

4. If £ is the "NNL rebuild" event then: 

(a) remove all events from calendar. 

(b) update the system to current time. 

(c) evaluate tr. 

(d) check for overlaps between parallelepipeds and build the NNL accordingly. 

(e) predict all cell-crossings and schedule them. 

(f) predict all collisions between ellipsoids and schedule them. 

(g) schedule next NNL rebuild. 

(h) schedule all other remaining events removed from calendar (output of summary, 
etc.). 

5. If £■ is a collision between particles A and B then: 



(a) change angular and center-of-mass velocities of A and B according to Eqs. (40) 

(b) remove from calendar all events (collisions, cell-crossings) in which A and B are 
involved. 

(c) evaluate new collision times Ia and of A and B with their parallelepipeds, 
and schedule a new "NNL rebuild" event if Ia or ts are less than tr. 

(d) predict and schedule the two cell crossings events for A and B. 

(e) predict and schedule all possible collisions for A and B using the NNL of A and 
B. 

6. If £ is a cell crossing of a certain rigid body A: 

(a) update linked lists accordingly. 

(b) remove from calendar all events (collisions, cell-crossings) in which A is involved. 

(c) predict and schedule all possible collisions for A using the NNL of A. 

7. Go back to step 2. 

Note that linked lists are used just to rebuild the NNL, that is given a parallelepiped 
corresponding to an ellipsoid A, one searches for overlapping parallelepipeds among all 
the ones that are in the same cell as A or in one of the 26 neighbors cells, assuming that 
the cells side is greater than the length of the diagonal of the parallelepipeds. To rebuild 
NNLs Donev et al. [26^ use, in addition to LL, also a collection of small spheres, called 
hounding sphere complex and in their implementation this method ensures that the cost 
of building the NNLs can be controlled increasing the elongation. In our implementation 
of NNL, that relies on the use of bounding parallelepipeds, the average collision time 
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(see discussion in Section 4.3) is quite independent from elongation for Xq up to 10. 
This is mainly due to the fact that the time needed to check overlaps between bounding 
parallelepipeds is negligible. 

In the above scheme NNL update is a global event (complete update) , in f2H| a method 
is described that allows to update only NNLs of affected particles after a collision (partial 
update) . In their implementation NNL are built using HEs instead of parallelepipeds and 
partial update method outperforms the complete update method (see ^6J). Differently 
in our implementation the time needed to update NNLs at moderate and high densities, 
which is what we are interested in the most, is negligible (see next Section) and the NNLs 
partial update method cannot offer any significant speedup. 

4. Performance results. 

To analyze the performance of the algorithm, monodisperse uniaxial HEs in the 
isotropic liquid phase have been simulated, characterized by the aspect ratio Xq — a/b 
(where a is the length of the revolution axis, b of the other ones) and by the packing 
fraction (f) = irXob'^p/G (where p = N/V is the number density). N = 256 particles 
have been simulated at various volumes V and elongations Xq in a cubic box with pe- 
riodic boundary conditions. The length of the smallest semi-axis is chosen as unit of 
distance and the mass m of the particle as unit of mass (m =1). Temperature is one in 
units of the Boltzmann constant ks', the corresponding unit of time is -y/ mP /{ksT). A 
spherically symmetric moment of inertia, i.e. Ix = ly = Iz = 2mr^ /b is chosen, where 
r = min(a, 6) /2 is the radius of the sphere inscribed in the ellipsoid. Notice that the 
value of the moment of inertia along the symmetry axis is not relevant for the present 
system, since angular velocity around the symmetry axis is zero. Although ellipsoids of 
revolution behave as symmetric tops, the colliding surfaces will be treated as perfectly 



smooth (see Eq. (40)) and therefore the component of angular velocity along the sym- 
metry axis will be conserved in each collision. In the following we indicate as "reduced 
time" the simulation time expressed in reduced units, while "CPU time" means the real 
time spent by the CPU for calculations, expressed in seconds. To test the algorithm 
speed, we use the CPU time per ellipsoid collision, labelled with Tc and defined as: 

rc=^ (69) 

I^coll 

where Ttot is the (real) time needed to perform Ncoii collisions during a simulation. 

Predicting a collision for a pair of colliding HEs (labelled by A and B) requires the 
CPU time McoU- 

AUn = X^<5if + X^(5tf (70) 

where N^^ indicates the number of collisions that have to be predicted for ellipsoid X 
{X = A,B) and with 6tf the CPU time requested to predict one collision. Indeed, after 
a collision between two HEs, all the possible future events involving the two colliding 
HEs must be predicted. If the quantity "average prediction time" is defined as follows: 

I 
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Figure 7: Time per collision Tc in the prolate case {Xq = 0.50) using only LLs (green diamonds) and 
also NNLs (orange circles) and in the oblate case (Xq = 2.00) using only LLs (red triangles) and also 
NNLs (blue squares). Dot-dashed lines are guides to the eyes, showing the <p and !/</) behavior for the 
NNLs case, while dashed lines are guide to the eyes showing a linear behaviour in <p for the LLs case 
(see discussion in Section |4.2| . NNLs show much less density dependence than LLs. Note that in all 
simulations the number of LL cells into which the box is partitioned is not kept fixed. As a consequence 
Tc vs <j> shows a simple linear dependence up to </> 0.55, where, due to increasing density, the number 
of LL cells decreases by one unit breaking the linear behavior of Tc(<ji). The inset is a zoom of Tc vs (f> 
employing also NNLs. 



where again X = A,B. Eq. (70 1 can be rewritten as 



Atco« = N^c{St)A + N^c{St)B (72) 

On passing note that, using LL and/or NNL, N^^ and iV^^ are 0(1) varying the total 
number N of HEs in the simulation. Such a number will be generally greater for LL 
than for NNL, but we will discuss this point more accurately later. Also note that {6t)A 
and {St)B depend mostly on the average number of Newton-Raphson steps performed 
to predict the collision time. In addition, during the simulation NNL and LL have to 
be updated and this will cost an extra time At^ per HE collision, so that the time per 
elHpsoid collision Tc is: 

T, = At^M + Atu (73) 

The contribution of the LLs update to At„ compared to the NNLs update contribu- 
tion is always negligible, meaning that the main contribution comes from the update of 
the NNLs, because it is more time consuming. Indeed, NNLs update requires the evalu- 
ation of the escape time of each HE from its bounding box and this is time consuming. 
Moreover, at high densities, At^ is usually negligible with respect to Atcoii, because the 
diffusivity of the particles is low at high densities and the average escape time of HE 
from their bounding boxes becomes quite long. All simulations are performed on a Intel 
Xeon CPU 3.06 GHz (codenamed "Prestonia") with a L2 cache size of 512 KB and a 533 
Mhz FSB. 



24 



4.I. Generation of the initial configurations 

To create the starting configuration at a desired 4>, an extremely diluted a-FCC crys- 
tal has been melted [58i ; afterwards, the particles have been grown independently up to 
the desired packing fraction (quench in at fixed N, Xq). The details of the growth algo- 
rithm will be illustrated in a future publication. To generate history-independent config- 
urations, tests have been performed on equilibrium configurations. To test the equilibra- 
tion of the sample, the decay of the self correlation function Fseif{q, t) — ^ e^' 
and of the orientational correlation function, C2(t) = P2{cos9{t)) - where P2{x) = {3x^ — 
1) /2 and 6(t) is the angle between the symmetry axis at time t and at time zero, have 
been studied. 



4-2. Dependence of the algorithm speed on density. 

Elongations Xq = 0.50 and 2.00 have been considered, fixing = 10^^ and the 
thickness of the neighbour shell Aat^vl = 0.15. 

Using only LLs both in the prolate and in the oblate case (see Fig. |7]) ), Tc increases 
going from lower to higher densities. In this case At^ is negligible and the average time 
to predict a collision (St) is roughly proportional to the number of HEs inside each LL 
cell, that in turn is proportional to the volume fraction (f> of the system. 

When using both LLs and NNLs the scenario is quite different. In Fig. [7] results 
from simulations using NNLs are shown, and it is apparent that there are two different 
regimes of Tc(0) at low and high volume fractions (f>. The interpretation of this observed 



behavior is quite straightforward. Indeed, considering Eq. (73 1, below w 0.5 there are 
few HEs within each NNL, i.e. N^^ « and N^^ « and the time per collision is 
dominated by Atu, that is by the time per collision needed to update the NNLs. In this 
case Atu is proportional to the number of average evaluations of the contact time of a 
HE with its bounding box tr, i.e. if I oc l/(f) is roughly the mean free path between two 
HEs collisions: 

Atu oc 1/Annl oc - (74) 



In Fig. [Tjthe ip~-^ dependence is evident at small (j). 

In contrast, at high densities Atu is negligible and the CPU time needed to predict a 
collision is dominated by AtcoU, that in turn depends linearly on the number of neighbors 
within the bounding box of a certain HE, i.e. 

47r 

Tcoc Naoc —(I) (a + ANNL){b + Annl){c + Annl) - 1 (75) 

At high volume fractions this equation simplifies to Tc oc and again Fig. [7] confirms this 
dependence on at high volume fractions. 

Finally, comparing NNLs and LLs results in Fig. [7) it is apparent that using NNLs 
the time per HE collision is significantly smaller than in the case where only LLs are 
used. Moreover, it is far less density-dependent. This stems from the fact that in Eq. 
N^^ and N^^ using NNLs are much smaller than the same quantities using LLs (see 
discussion in the next section). 
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Figure 8: Dependence of average coUisione time Tc on elongation Xq. (a) at <f> = 0.30,0.40 and 
0.50 using only LLs. In this case the average collision time is linked to the number of particles in the 
first neighbour shell and grows with elongation. The dashed lines are guides to the eyes showing the 
1/Xo and Xq behaviors (see text), (b) Comparison of performance results with and without the use of 
NNLs plotting ^/r^^-^ versus the aspect ratio Xq for 3 different volume fractions </) = 0.30, 0.40, 0.50. 
^NNL jg ^]jg collision time employing NNLs and t^^ is the collision time using only LLs. The use of 
LLs severely degrades the performance of the algorithm for elongated HEs. 



4-3. Dependence of algorithm speed on elongation 

In the following we will show results for simulations performed at volume fractions 
(j) — 0.30,0.40,0.50. First, the case where only LLs are employed is considered. Besides 
the expected decrease of Tc upon increasing the density, a marked increase of Tc upon in- 
creasing/decreasing the elongation with respect to the hard-sphere case Xq = 1 is clearly 
apparent in Fig. |8]^a). The explanation for such behavior is straightforward: as noticed 
in [5^, the number of HEs in a LL cell increases as Xq (I/Xq) for the prolate (oblate) 
case, and in turns Tc increases as Xq {1/Xq) at big (small) elongations. In contrast, 
using NNLs Tc is quite Xo-independent, increasing significantly the performance at big 
and small elongations with respect to LLs (see Fig. |8jb)). In this case the number of 
neighbors is independent from on changing the elongation and the computational effi- 
ciency depends only on NR method efficiency, that is quite Xo-indipendent. On passing 
we note that at very high elongation ( Xq ^ 5), the NR method requires more steps to 
converge and the exact number of steps depends on the initial guess. Nevertheless, we 
checked that, at least up to Xq = 10, the algorithm's performance does not change (i.e. 
Tc remains Xp-independent) if NNLs are used. 



4-4- Optimizing the parameter 



The main parameter entering in the collision prediction algorithm is e^j (sec. 2.4.2) 



and it is important to establish the optimal value in terms of efficiency. As already 
discussed in 2.4.2 too large a value for this parameter may result in a failure of EDMD 
due to overlaps of HEs. Hence, the best choice for ed is the largest value not generating 
HEs overlaps. Fig. |9] shows that Tc monotonically decreases as is decreased, and that 
^dependence is rather weak. In practice, a value like = lO""* — 10~^ is usually a good, 
reliable and safe choice for most situations. 
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Figure 9: Dependence of Tc on (using NNLs) for two different elongations {Xq = 0.5, 2.0) at various 
values of ranging from 0.30 to 0.55. Arrows indicate decreasing volume fractions. 



4-5. Optimizing neighbour lists 

Changing the parameter A^nl, i.e. changing the size of the bounding boxes, compu- 
tational times vary non-monotonically. If Annl tends to 0, the number of neighbors N^^, 
Np^ in Eq.(72l tends to accordingly, but NNLs updates tend to be infinitely frequent, 
because escape time of HE from their bounding boxes tends to 0. Therefore diverges 
to infinity, if Ajvatl 0. If A^vatl ~^ oo (assuming one has a system with infinite par- 
ticles), escape times tends to and time intervals between two successive NNLs rebuilds 
diverge, but N^^ oo, iV^^ — > oo in Eq.(72l, so that again Tc ^ oo. As a result, there 
must be a minimum of Tc as a function of A^r jvl • The value of A nnl , 
Tc is the optimal value for this parameter and it will be labeled A^^^. 



which minimizes 
Fig. 



10 shows for 



two volume fractions that Tc(Ajvjvl) exhibits a clear minimum. A^^^ depends much 
more on the volume fraction (p than on the elongation Xq and this is clear from Fig. [TT] 

for different Xq. The latter result is due to the 



where A^jy^ is plotted as a function of 



pc 1 



A^^, strongly depends on (j), while both At 



coll 



fact that in Eq. (72 1 Atcoii, through 
and Atu weakly depend on X^. 

Fig. 11 provides a simple estimate of A^^^, that ranges from 0.1 to 0.8 if 0.5 < 
Xq < 2.0 and 0.25 < < 0.55, upon changing volume fraction and elongation. 



5. Simulating Superquadrics. 

In the Sections|3]and[4]we have shown in detail how to apply the algorithm illustrated 
in Sec. [2] to simulate HEs and its performance in this particular case. The algorithm 
proposed by Donev et al. |25} I26j has been generalized to arbitrary convex shapes 
PSI [301 E], analogously we consider here the simulation of SQs, that are a possible 
generalization of HEs studied in the previous section, skipping anyway some details that 
will be supplied in a future publication. 
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Figure 10: The average collision time depends non-monotonically upon A^vjvLi showing a clear minimum 
for a certain AJ^^^ value (see text). In this figure such a dependence is shown for 4 particular state 
points (Xo = 0.5 a.t (j) = 0.45, 0.57 and Xq = 2.0 at = 0.45, 0.57). 



5.1. Geometry and Simulation Details 



The surface of an HE has been defined through Eq. (42 1, that with a suitable choice 



of the reference system axes (principal axes) takes the form: 

c 



^ a / 



ZX2 



1 



(76) 



A straightforward generalization of the latter equation leads to the definition of particles, 
called SQs, whose surface is the following: 
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1 = 



(77) 



where the parameters n,m,p are real numbers and a, b, c are the SQ semi-axes. A 
monodisperse system of N = 512 SQs has been simulated with m = n = 2 and with two 
equal semi-axes, i.e. b — c. 

Such SQs can be characterized by the elongation Xq = a/b and by the parameter p, 
that determines the sharpness of the edges (see Fig. 12 1. SQs of elongation Xq < 1 are 
called "oblate" , while SQs of elongation Xq > I are called "prolate" , as for the HEs. 
Elongations Xq = 0.5 and Xq = 2.0 have been studied. Nicely for Xq = 0.5 and p = 8 
the SQ has a pillow-like shape, while for Xq = 2.0 and p = 8 it has a boat-like shape, as 
shown in Fig. 12 The system of prolate and oblate SQs has been simulated in a cubic 
box of volume V with periodic boundary conditions at the volume fraction (j) = 0.256 
and at various values of p. The length of the smallest semi-axes is chosen as 0.8 times 
the unit of distance, the mass of the SQ as unit of mass (to = 1) and the moment of 
inertia is spherically symmetric, i.e. Ix ^ ly = Iz = 1, as in the case of HEs. 

The code for simulating SQs can be straightforwardly derived from the HE code, 
indeed for implementing the EDMD of SQs it suffices to: 



Evaluate the Jacobians in Eqs. ^ ( or Eqs. (23l ) and (31 1 using Eq. (77) 
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Figure 11: In this figure the dependance of A^^y^ on for Xq = 0.5, 2.0 (left panel) and the dependance 



of A^jY^ on Xq for (p = 0.40, 0.55 (right panel) is shown. A^ 
(left panel) but it has a weak dependence on Xg (right panel). 



decreases significantly with increasing 



Calculate the Jacobians pertaining to the collision of a SQ with a plane to make 



use of NNL (see Sec. 3.3 1. 



• Adapt the steepest descent method described in Sec. 2.3.2 

• Adapt the guess of the distance shown in Sec. |3.2[ 



5.2. Performance results 

Fig. 13 shows the average collision time Tc for the system of SQs using either LLs or 
NNLs. With regard to the dependence of Tc on (f>, Xq, td and A^nl We do not expect 
a significantly different behavior from what has been observed for HEs. Hence here we 
focus on the dependence of the average collision time on the parameter p. 

For the LLs case p ranges from 2 to 8, while enabling NNLs it ranges from 2 to 4 only. 
In the case of LLs for values oi p > 8 the NR for calculating the distance between SQs 
starts having convergence problems, while in the case of NNLs the NR for calculating 
the distance between a SQ and a plane starts having convergence issues. In the present 
implementation of the algorithm for locating the collision time between SQs and between 
a SQ and a plane the greater the p the slower will be the NR method to reach convergence 
(up to a maximum value above which the NR method starts failing) and the greater will 



be the number of steps needed to locate the contact point (see Sec. 2.4.21. 

About the performance comparison between HEs and SQs from Fig. |13| it turns out 
that the SQs simulations (see p = 2, i.e. HEs) are at least 5-6 times slower than the HEs 
ones and this is mainly due to the more time consuming calculation of the Jacobians 
in Eqs (23) and ( pTj ). Finally, as discussed above it is apparent from Fig. 
increase of Tc at high p in all cases. 
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6. Possible applications of the present algorithm 

The algorithm illustrated in this paper has been already applied in several fields, like 
biophysics, soft-matter and condensed-matter physics. First of all, this new algorithm 
was employed to investigate the dynamics of HEs [43,, finding, for the first time, evidence 
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Figure 12; Shape of the superquadrics used in the simulations. Top: on the left a prolate SQ with p = 8 
and Xq = 0.5 is represented in 3D. On the right the plots of the contour resulting from the intersection 
of the SQ (for several p ranging from 2 to 8) with the plane x = are shown. Bottom: on the left an 
oblate SQ with p = 8 and Xq = 2.0 is represented in 3D. On the right the plots of the contour resulting 
from the intersection of the SQ (for several p ranging from 2 to 8) with the plane y = are shown. 
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Figure 13: Average collision time Tc for simulations of N = 512 SQs at c/) = 0.256 for elongations 
Xo = 0.5 and Xg = 2.0 using only LL (left panel) and also using NNL (right panel). 
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of a pre-nematic order-driven glass transition. The fact that this algorithm can simulate 
hard objects of convex arbitrary shape offers new intriguing possibilities for investigating 
the changes in the phase diagram of HEs on changing particles shape. For example, 
concerning the pre-nematic order driven glass transition observed in [43] , it is stimulating 
to think about the possibility to completely inhibit the nematic transition with a different 
choice of the HRB shape and/or introducing shape polydispersity. Work in this direction 
is under way. Changing the HRB shape may result in a phase diagram, which is richer 
than the one of simple HEs, since new phases (cubatic, smectic, ecc.) may appear [5]. 
Moreover, changing the shape and introducing shape polydispersity may have effects on 
packing and on mixing/demixing of hard objects. 

This algorithm has also been extended j44] to combinations of hard and localized 
attractive interactions, modeled as site-site square well potentials. In this way, the hard- 
body may be decorated with an arbitrary number of spherical patches (sticky spots) 
arranged in fixed site locations. The algorithm for the prediction of the collision time 
between "sticky spots" is an adaptation of the algorithm illustrated in the present paper 
for HEs (some details can be found in [44 but this topic will be discussed extensively in 
a forthcoming publication) . 

Spheres decorated with sticky spots have been used to simulate primitive models of 
water and silica "^S] and to study the kinetics of self-assembly of an idealized fluid 
model |48j . Hard eUipsoids with sticky spots have been employed to study a model of a 
chemical gel [46] . 

HEs decorated with attractive interacting sites may also be used to model complex 
hard molecules, retaining some degree of flexibility. In this approach, the loss in the 
details of the system is counterbalanced by its extreme flexibility and by the possibility 
of investigating time scales which are not accessible by standard methods (like "ab-initio" 
calculations or full-atom molecular dynamics simulations). There are several directions in 
which this methodology could be applied. A very promising fleld of interest is represented 
by biophysical systems. A particular problem that has received much attention over the 
last years is the IgG antibody-antigen interaction. An immunoglobulin (or antibody) is a 
"Y" shaped protein that is ubiquitous in most vertebrates. Its role is to bind to viruses or 
bacteria to facilitate their neutralization process. It is extremely interesting that a single 
object has the ability to bind to biological targets whose size can vary from much smaller 
to much larger than its own. The introduction of square well interactions |44j between 
hard-body particles can be used to generalize the Go model [59 and better account 
for steric effects between different proteins. Previously, primitive model of proteins had 
been developed to study biophysical problems. In the 4-beads model, described in [SU] . 
each amino acid is represented by a maximum of four beads. Three beads correspond 
to the amide N, the a-carbon C, the carbonyl C groups. The fourth bead models the 
amino-acid side-chain group of atoms, and it is placed at the centre of the nominal 
atom. Exploiting the new algorithm proposed in the present paper, each amino-acid 
with its interacting (active) sites can be modeled as a set of spherical spots that form 
one unique rigid body. In biological simulations it is worth noting that the solvent can 
be also introduced explicitly in several ways: 

• using a primitive model, e.g. the primitive model of water used in [44] to take into 
account the directionality of the hydrogen bonding. 

• using a simplified solvent of hard spherical particles (this technique has been used 
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for studying the interaction of IgG and antigenes) 

• adapting the method for performing brownian dynamics of square-well particles 
discussed in [SH] to the present algorithm. 

Finally, this algorithm may find applications in the field of computer science: ani- 
mations and virtual reality may benefit from the physical accuracy and efficiency of the 
present method. In general, virtual motions of objects in 3D space require an accurate 
prediction of collisions and are fundamental for robotics, computer graphics, 3D com- 
puter games [61 and CAD applications. For example, because of its flexibility in shape, 
ellipsoids are often chosen as bounding volumes for robotic arms in collision detection 



7. Conclusions 

In this paper a new efficient method for performing event-driven molecular dynamics 
simulations of non-spherical HRBs has been proposed. This new method is based on the 
traditional NR method for solving a set of non-linear equations. In particular, geometrical 
distance and contact point and time between two moving HRBs can be calculated through 
the NR method in a very efficient way. The method has been tested against a well 
established and studied system, the hard ellipsoids of revolution and also against a less 
common system like the SQs. Anyhow HEs and SQs are just particular cases for the 
present algorithm, since more general convex hard bodies can be simulated. Furthermore 
also non-convex hard bodies can be in principle simulated by the present algorithm, if 
all the possible solutions of Eqs. ^ (or Eqs. (23)) can be found, because in this case 



the actual distance can be simply obtained by Eq. (11 1. 

In the specific case of uniaxial HEs the time per collision achieved at moderate and 
high densities for elongations Xq — 0.5, 2.0 is around 2.5 ms, a value which is comparable 
to the performance of the method described in |26j . For comparison, in our implementa- 
tion the time per sphere collision is about 0.25 ms, i.e. simulating hard spheres is about 
an order of magnitude faster than simulating HEs, even for nearly spherical HEs (the 
same observation has been carried out in [26 ). 



To solve the set of equations to evaluate the contact point and time (see Eqs. ( 30 )) by 



the NR method, a good initial guess ("bracketing") of the solution has to be provided. It 



has been shown (see 2.4.21 that an initial guess for Eqs. (30) can be found by evaluating 
the geometrical distance between the two colliding HRBs and by making use of a simple 
overestimate of the rate of variation of the distance (dmax)- It is worth stressing that 
using such an algorithm for bracketing the contact point and time, "grazing collisions" do 
not constitute a problem within machine accuracy. If Cd ~ 10~* all collisions are correctly 
predicted and this choice of ed does not depend on volume fraction or aspect ratio of 
simulated HEs. Also a new method to implement NNL based on oriented bounding 
boxes has been presented. This new NNL method is fast and flexible enough to be easily 
extended to more complex shapes than simple HEs and SQs. 

The present algorithm can also be generalized, without losing numerical efficiency, 
to the case of attractive interactions modeled via discontinuous step-wise potentials, 
including the interesting case of site-site square well interaction, providing the possibility 
of modeling highly directional interactions between the rigid bodies. This possibility can 
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be further refined by transforming site-site interactions into a p(^rniancnt link between 
objects, by simply changing the depth of square well interactions. In this way, the objects 
can be linked into flexible structures. Such flexibility couples well with the novel NNL 
implementation described. 
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